home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / DBRENT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  3KB  |  112 lines
  1. FUNCTION dbrent(ax,bx,cx,tol: real; VAR xmin: real): real;
  2. (* Programs using routine DBRENT must define the external
  3. functions func(x:real):real and dfunc(x:real):real which are,
  4. respectively, the function whose minimum is to be found,
  5. and its derivative. *)
  6. LABEL 1,2,3;
  7. CONST
  8.    itmax=100;
  9.    zeps=1.0e-10;
  10. VAR
  11.    a,b,d,d1,d2: real;
  12.    du,dv,dw,dx: real;
  13.    e,fu,fv,fw,fx: real;
  14.    iter: integer;
  15.    olde,tol1,tol2: real;
  16.    u,u1,u2,v,w,x,xm: real;
  17.    ok1,ok2: boolean;
  18. FUNCTION sign(a,b: real): real;
  19.    BEGIN
  20.       IF (b > 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
  21.    END;
  22. BEGIN
  23.    IF ax < cx THEN a := ax ELSE a := cx;
  24.    IF ax > cx THEN b := ax ELSE b := cx;
  25.    v := bx;
  26.    w := v;
  27.    x := v;
  28.    e := 0.0;
  29.    fx := func(x);
  30.    fv := fx;
  31.    fw := fx;
  32.    dx := dfunc(x);
  33.    dv := dx;
  34.    dw := dx;
  35.    FOR iter := 1 TO itmax DO BEGIN
  36.       xm := 0.5*(a+b);
  37.       tol1 := tol*abs(x)+zeps;
  38.       tol2 := 2.0*tol1;
  39.       IF (abs(x-xm) <= (tol2-0.5*(b-a))) THEN GOTO 3;
  40.       IF (abs(e) > tol1) THEN BEGIN
  41.          d1 := 2.0*(b-a);
  42.          d2 := d1;
  43.          IF (dw <> dx) THEN  d1 := (w-x)*dx/(dx-dw);
  44.          IF (dv <> dx) THEN  d2 := (v-x)*dx/(dx-dv);
  45.          u1 := x+d1;
  46.          u2 := x+d2;
  47.          ok1 := ((a-u1)*(u1-b) > 0.0) AND (dx*d1 <= 0.0);
  48.          ok2 := ((a-u2)*(u2-b) > 0.0) AND (dx*d2 <= 0.0);
  49.          olde := e;
  50.          e := d;
  51.          IF (NOT (ok1 OR ok2)) THEN GOTO 1
  52.          ELSE IF (ok1 AND ok2) THEN BEGIN
  53.             IF (abs(d1) < abs(d2)) THEN BEGIN
  54.                d := d1
  55.             END ELSE BEGIN
  56.                d := d2
  57.             END
  58.          END ELSE IF (ok1) THEN BEGIN
  59.             d := d1
  60.          END ELSE BEGIN
  61.             d := d2
  62.          END;
  63.          IF (abs(d) > abs(0.5*olde)) THEN GOTO 1;
  64.          u := x+d;
  65.          IF (((u-a) < tol2) OR ((b-u) < tol2)) THEN BEGIN
  66.             d := sign(tol1,xm-x)
  67.          END;
  68.          GOTO 2
  69.       END;
  70. 1:      IF (dx >= 0.0) THEN e := a-x ELSE e := b-x;
  71.       d := 0.5*e;
  72. 2:      IF (abs(d) >= tol1) THEN BEGIN
  73.          u := x+d;
  74.          fu := func(u)
  75.       END ELSE BEGIN
  76.          u := x+sign(tol1,d);
  77.          fu := func(u);
  78.          IF (fu > fx) THEN GOTO 3
  79.       END;
  80.       du := dfunc(u);
  81.       IF (fu <= fx) THEN BEGIN
  82.          IF (u >= x) THEN a := x ELSE b := x;
  83.          v := w;
  84.          fv := fw;
  85.          dv := dw;
  86.          w := x;
  87.          fw := fx;
  88.          dw := dx;
  89.          x := u;
  90.          fx := fu;
  91.          dx := du
  92.       END ELSE BEGIN
  93.          IF (u < x) THEN a := u ELSE b := u;
  94.          IF ((fu <= fw) OR (w = x)) THEN BEGIN
  95.             v := w;
  96.             fv := fw;
  97.             dv := dw;
  98.             w := u;
  99.             fw := fu;
  100.             dw := du
  101.          END ELSE IF ((fu < fv) OR (v = x) OR (v = w)) THEN BEGIN
  102.             v := u;
  103.             fv := fu;
  104.             dv := du
  105.          END
  106.       END
  107.    END;
  108.    writeln('pause in routine DBRENT - too many iterations');
  109. 3:   xmin := x;
  110.    dbrent := fx
  111. END;
  112.